Bound states and E 8 symmetry effects in perturbed quantum Ising chains 
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In a recent experiment on C0ND2O6, Coldea et al. [IJ found for the first time experimental 
evidence of the exceptional Lie algebra Eg. The emergence of this symmetry was theoretically 
predicted long ago for the transverse quantum Ising chain in the presence of a weak longitudinal field. 
We consider an accurate microscopic model of C0ND2O6 incorporating additional couplings and 
calculate numerically the dynamical structure function using a recently developed matrix-product- 
state method. The excitation spectra show bound states characteristic of the weakly broken Eg 
symmetry. We compare the observed bound state signatures in this model to those found in the 
transverse Ising chain in a longitudinal field and to experimental data. 

PACS numbers: 75.10.Pq, 75.40.Mg, 75.78.Fg 



The one-dimensional (ID) quantum Ising model in 
transverse and longitudinal fields is one of the most stud- 
ied theoretical models in condensed matter physics. It is 
a relatively simple model that contains very rich physics; 
for example, it contains a quantum critical point (QCP) 
at zero longitudinal field related to the 2D classical Ising 
model. A remarkable fact is that the integrability present 
at the critical point remains under addition of a longitu- 
dinal field as a mass-generating perturbation. Zamolod- 
chikov conjectured in 1989 an S-matrix describing eight 
emergent particles whose mass ratios are connected to the 
roots of the Lie algebra Eg [H, Q . Recently, Coldea et al. 
performed neutron scattering experiments on CoNb2 06 
(cobalt niobate), a material that to a good approxima- 
tion can be described by a quantum Ising chain. At low 
temperatures and in the presence of a strong external 
transverse magnetic field which tunes the system to near 
criticality, the observed spectrum shows characteristic ex- 
citations of the Eg symmetry [l[ . 

However, a serious problem in comparing theory and 
experiment is that the real material has additional cou- 
plings that strictly speaking invalidate the exact solution, 
and until recently it was impractical to extend the the- 
ory non-perturbatively to include these couplings. In this 
Letter, we study a theoretical model for CoNb206 which 
includes in addition to the Ising interaction other inter- 
actions arising from the lattice structure and the weak 
coupling between the chains. Using this model, we cal- 
culate the dynamical spectral function and compare the 
results to the observed spectra. Close to the QCP, the 
model retains features expected from the quantum Ising 
model, in particular the characteristic particles of the Eg 
symmetry. 

We begin by deriving the theoretical model used to 
describe the low-energy physics of CoNb20g. The spin 
lattice structure consists of chains of easy axis spins, re- 
alizing a two level system, on the Co 2+ ions coupled by a 
ferromagnetic Ising interaction along the chain direction, 
see Fig.QJ^. We thus start from the quantum Ising chain, 



described by the Hamiltonian 

H = — J S z t S z l+1 — h x S x 



(1) 



where J > favors a ferromagnetic state (| fj ... f) 
or I 11 ... I)). When the transverse field is increased 
past the QCP \h x \ = J/2, the system undergoes a phase 
transition into a paramagnetic state | — >^ ... — >). This 
model is exactly solvable using a Jordan- Wigner transfor- 
mation which transforms the spins into non-interacting 
fermions 

The lowest lying excitation energy is similar on both 
sides of the QCP due to the self-duality of the model and 
goes to zero, that is, the gap closes, at the QCP. How- 
ever, the double degeneracy of the ferromagnetic ground 
state leads to a fractionalization of the experimental ex- 
citation, a spin flip, into two freely moving domain walls 
or kinks. We now take into account terms which re- 
sult from the three-dimensional (3D) lattice structure of 
CoNb2 0g. A recent theoretical study by Lee et al. inves- 
tigates a three-dimensional model of CoNb206 They 
show that the plane perpendicular to the chain, a weakly 
coupled triangular lattice, see Fig. [TJA., has ferrimagnctic 
order to transverse field strengths well passed h x . The 
interchain couplings in a 3D magnetic ordered material 
at low temperature can be well approximated by a chain 
in a local effective longitudinal field h z = y\ J$(S Z ) with 
the sum over all nearest interchain bonds |6|. This field 
favors the ferromagnetic phase, breaks its two-fold sym- 
metry and moves the system away from the QCP. It also 
splits up the continuum into bound states by confining 
the kinks. At low transverse field and small bound state 
momentum, this can be described by a one-dimensional 
Schrodinger equation with a linear confining potential 
with the energy levels given by the negative zeros of the 
Airy function, see Fig. [TJ3 [7[. This solution has later 
been extended to all possible bound state momenta [§). 
Close to the QCP (h x = h x , \h z \ < \h x \), the eight mas- 
sive particles described by the Es symmetry can be seen 
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either as asymptotic states or as bound states of a pair 
of particles of this theory 0, Q . 

Although CoNb 2 06 to a good approximation can be 
described by a quantum Ising chain, a realistic model 
must contain more interactions Q. It has a strong easy 
axis character, but a weak XX part is still present. The 
chains have a zig-zag structure, making the next-nearest 
neighbor (nnn) interaction important as well, see Fig.[JJ\. 
The measured Ising exchange energy J is unusually low, 
likely due to a competition from an antiferromagnetic 
nnn interaction. Taking into account all terms, the re- 
sulting Hamiltonian reads 

H = -J'^2S*S* +1 - h x ^S x - h'^^S* (2) 

n n n 



the truncation effects on local expectation values (< 
for all simulations presented in this paper). In addition 
we also checked the dependence of the measured observ- 
ables on the matrix dimension x P er site and the time 
steps At, settling for \ = 45 and At = 0.04 meV -1 for 
the simulations presented in this paper. 



Schematic structure of Co ions 
. in CoNb 2 O e 
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The numerical values for the coupling constants to de- 
scribe CoNb206 are obtained by matching the experi- 
mental neutron scattering intensity at zero applied trans- 
verse field with our numerical calculations. We compare 
the dynamical structure function S v (k,u)), the Fourier 
transform of the dynamic two-point correlations 



Cy(n,t) = {i> \SZ(t)S%(0)\iJo). 



(3) 



For the numerical calculations, 



we use the time evolv- 



ing block decimation (TEBD) [id . 1 1 11 ] method which pro- 
vides an efficient method to perform a time evolution of 
quantum states in one-dimensional systems. The evolu- 
tion of a random state of an infinite chain in imaginary 
time is used to calculate the ground state \tpo) and an 
evolution in real time allows us to calculate the dynamic 
two-point correlations directly. The TEBD algorithm 
can be seen as a descendant of the density matrix renor- 
malization group [l2| method and is based on a matrix 
product state (MPS) representation [nl, [3] of the wave- 
functions. Algorithms of this type are efficient because 
they exploit the fact that the ground-state wave func- 
tions are only slightly entangled, especially away from 
criticality [15j. As the entanglement grows linearly as a 
function of time, the simulations of long time evolutions 
is numerically very difficult. To be able to simulate long 
enough times and thus to get sufficiently good energy 
resolution in the calculated spectral functions, we use a 
number of methods to accelerate the time evolution. We 
use linear predictions to extrapolate the d yna mical cor- 
relation functions to very long times [1(3, 17 1 and take 
advantage of the "light-cone" like spread of the entangle- 
ment by adding more sites to the chain as time increases. 
As the calculation of the correlation functions C v (n, t) is 
numerically very expensive, we calculate it only for cer- 
tain time steps and then interpolate its values. In order 
to estimate the errors of our simulations, we calculate the 
truncation error, i.e., the truncated weight of the wave 
function at a time step, which gives an upper bound for 
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FIG. 1: (Color online) (A) Ising spins on the Co 2+ ions are 
strongly coupled in ID along zig-zag chains. The Co 2+ ions 
are ordered in a weakly coupled triangular lattice in the plane 
perpendicular to the chain direction, with a, b and c orthogo- 
nal unit vectors. (B) Confinement of kinks: the potential en- 
ergy between kinks increase linearly, along the z coordinate 
in the c direction, as more interchain bonds turn energeti- 
cally unfavorable. The energy levels are given by the nega- 
tive zeros of the Airy function. (C-D) The full Hamiltonian 
describing CoNb2 0e at no external magnetic field. (C) The 
dynamical structure function. (D) The cross section of (C) 
at zero momentum showing the masses of the first five bound 
states. A comparison of these masses from our MPS calcula- 
tions (pluses) with the experimental results (circles) and the 
exact solution of the proposed first order phenomenological 
model (crosses). 



The numbers we use are J' = J + Jb = 2.43 meV, 
h x = 0.354 meV, h* = 0.035 meV, J p = 0.52 meV and 
Jb = 0.60 meV, see Fig. [Tp and compare it to Fig. 3 
in Ref. [l|. A cross-section with the bound state "masses" 
(i.e., the energies of the bound states at zero momentum) 
is presented in Fig. [Tp with the experimentally mea- 
sured masses and Rutkevich's exact solution of Coldea 
et at first order model for reference in Fig. [Tp [l|, [lj|. 
Note that our full Hamiltonian agrees to first order in per- 
turbation theory with the phenomenological model used 
there, although the coupling constant for that model is 
slightly larger than ours 
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FIG. 2: (Color online) Ising chain in a transverse and a 
longitudinal field. (A) The dynamical structure function at 
h x = h% and fe z = 0.035 meV. (B) The cross section of (A) at 
zero momentum. The five lowest bound states and two bound 
state pairs can be distinguished. (C) (h x = h x ) The mass of 
the four lowest bound states and an asymptotic expansion for 
mi from h z — > as a function of the longitudinal field. (D) 
(h x = h*) Relative mass of the lowest bound states compared 
to the analytical predicted values at h z — ¥ 0. (E) (h z = 0.035 
meV) The bound state masses as a function of the transverse 
field (h x ). The energy gap is smallest around h x — 1.10 meV, 
however the minimum in the higher bound states occurs for 
lower fields. (F) (h z = 0.035 meV) The relative masses for 
the bound states increases roughly linear as a function of the 
field (h x ) and passes the analytical calculated values at h x . 



Before investigating this model near the QCP, we start 
with a pure quantum Ising chain where the Es symmetry 
is expected to be present The dynamical structure 
function is calculated over the whole Brillouin zone (for 
various parameters J = 1.83 meV, h x , h z ), see Fig. [2j\ 
for an example and we focus on the cross section at zero 
momentum where some comparison with earlier work can 
be done. Fig. 03 is an example of these cross sections at 
h x = h x with h z = 0.035 meV. The lowest four bound 



states can be easily detected and one or two more can 
be distinguished. Bound state pairs mi + mi (overlap- 
ping with 7713) and mi + 777,2 have similar intensity to the 
nearby bound states, making both types simultaneously 
observable. They are created in a "spinon jet", where the 
two kinks (also known as spinons) in a bound state have 
been stretched far enough apart to make it energetically 
favorable to flip a spin between them to form two more 
kinks that each form an independent low energy bound 
state with one of the original kinks. The independent 
motion of these two bound pairs appears as a contin- 
uum in the dynamical structure function. This process 
is reminiscent of quark dynamics, where the quarks are 
confined and cannot be isolated singularly. Finding con- 
densed matter analogues of confinement effects known 
from high energy physics might help us to improve our 
understanding of underlying mechanisms; see e.g. Lake 
et al. [111. 

The weight of the continuum decreases with increasing 
longitudinal field; this is also the case for the weight of 
the higher bound states but to a lesser extent. The gap 
and the spacing between the bound states increase with 
increasing longitudinal field; sec Fig.[2jZI where data from 
more simulations are presented, together with an asymp- 
totic expansion from the exact analytical limit (h x = h x , 
h z — > 0) of the lowest bound state. The analytical ex- 
pression for the gap is mi « CJ/4{2h z /J) 8 / 15 , with 
C = 4.40490858/0.7833, showing good agreement with 
our results to high longitudinal field strengths [H, [2fj| . 
(The spin is rescaled Sf at (x) = 0-783(3}S z ont (x) from the 
continuum to the lattice model [U [22[.) The relative 
mass of these bound states related by the Eg symmetry 
at he, are presented in Fig. 0D, again with good agree- 
ment in the analytically exact limit, see Tab. U 

7712 /mi 7773 /mi 7714/mi Tils /mi me/mi my /mi ms/mi 
1.618 1.989 2.405 2.956 3.218 3.891 4.783 

TABLE I: Analytically predicted mass ratios from Ref. 0. 
These numbers result from evaluation of simple trigonometric 
expressions (e.g., mi/mi = 2cos7r/5) that arise as eigenval- 
ues of a matrix constructed from roots of the Lie algebra Es. 



The deviation of the asymptotic expansion is slightly 
larger for higher bound states. However, the deviation for 
large longitudinal fields (h z < h x ) is fairly small, indicat- 
ing influence of criticality up to very strong longitudinal 
fields. Note also that the uncertainty of our results in- 
crease with decreasing longitudinal field strength when 
the bound state masses move closer, due to our fixed 
energy resolution. For future reference we also present 
results as a function of the transverse field at longitudi- 
nal field h z = 0.035 meV present in C0NLJ2O6 around 
h x = h x , see Fig. Efl. Good agreement for the energies 
of the bound states are obtained with earlier numeri- 
cal work using the Truncated Free Fermion Space Ap- 
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proach, cf. Fig. 5 of Rcf. |23j. Stronger longitudinal field 
will increase the minimum gap and move it to stronger 
transverse fields, but the gap increases slower away from 
its minimum value. Also note that the minimum for 
higher bound state masses occurs for a lower transverse 
held. The relative masses increase linearly around h*, 
see Fig. [2F, with a steeper slope for higher bound state 
masses and lower longitudinal fields. 
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FIG. 3: (Color online) The full Hamiltonian describing 
C0ND2O6. (A) Magnetization comparison between weakly 
coupled chains and chains in different constant longitudinal 
field. (B) The cross section of the dynamical structure func- 
tion for h% at zero momentum showing the masses of the 
first five bound states and two bound state pairs. (C) The 
bound state masses as a function of h x . The minimum gap is 
above h x c and the bound state mass minimum decreases with 
increasing mass. (D) The ratio of the bound state masses 
varies linearly around and goes through the analytically 
calculated values at h%. 



the one for the quantum Ising chain, with fairly unaltered 
spectral weights. The relative weight of the bound state 
continuum is still largest around h*, making this region 
even more interesting for experiments. A more careful 
analysis of the bound state masses, see Fig. reveals 
a small rescaling of both axes to around 90% of their 
previous values. This overall scaling does not affect the 
mass ratios, see Fig.|3p. Again they follow straight lines, 
see Fig. 03, and pass the analytical values at h*, exactly 
as they do for the quantum Ising chain, rather than ap- 
proaching it by bending as suggested by the extrapolation 
of the experimental data in Ref. [l|. Additional interac- 
tions irrelevant at low field and not treated here might 
explain the bending if it is confirmed by higher resolution 
data, but our results suggest higher resolution data will 
show that the mass ratios indeed go through the analyti- 
cal values at the critical field if the model used here (and 
previously [H, 0]) is a good one for CoNL^Og past h x . 

To conclude, we have investigated the effects of inte- 
grability near the Ising QCP and evaluated how far away 
the features extend and how robust they are to addi- 
tional interactions. We have shown that the bound state 
continuum should carry comparable spectral weight to 
the higher bound states. The microscopic ID model of 
CoNb 2 06 treated here is able to reproduce the experi- 
mental data far away from criticality well. When moved 
close to the QCP, the model still has the characteristics of 
the Eg symmetry, with the mass ratios following straight 
lines through the analytical values, even better than the 
extrapolated experimental data suggests. Future experi- 
ments with improved resolution on CoNb2 0g should de- 
tect higher bound state signatures to confirm the effects 
of integrability and the bound state continuum modeling 
confinement dynamics around the ID QCP. 
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OLE program and by the Knut and Alice Wallenberg 
foundation (J. K.). 



Finally we turn to the more accurate microscopic 
model of CoNb 2 06 Eq. with values of the coupling 
constants presented above. The QCP at zero longitu- 
dinal field for this model is moved to a slightly weaker 
field h x « 0.814 meV, see Fig. [3]A. from ground state 
simulations with TEBD, due to the addition of the ferro- 
magnetic XX-term. The longitudinal field strength from 
weakly coupled chains h z ((S z )) is to a good approxima- 
tion constant past h*, see Fig.[3]A At vanishing magneti- 
zation this is not true, but here the ID approximation of 
the 3D material is breaking down anyway. The dynamical 
structure function at h x , not presented here, shows a flat- 
tening of the kinetic bound state and a more prominent 
lowest bound state dispersion. The cross section at zero 
momentum, see Fig. 03, has the same characteristics as 
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